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We illustrate the scope of Time Dependent Density Functional Theory (TDDFT) for strongly 
correlated (lattice) models out of equilibrium. Using the exact many body time evolution, we 
reverse engineer the exact exchange correlation (xc) potential v xc for small Hubbard chains exposed 
to time-dependent fields. We introduce an adiabatic local density approximation (ALDA) to v xc 
for the ID Hubbard model and compare it to exact results, to gain insight about approximate xc 
potentials. Finally, we provide some remarks on the v-representability for the ID Hubbard model. 
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Density Functional Theory (DFT) [JJ enables accurate 
investigations of realistic systems of considerable com- 
plexity. However, strongly correlated systems (SCS) have 
until now remained elusive to DFT. And, so far, most 
theories of SCS focus on equilibrium or non-equilibrium 
steady state regimes, to understand the long time re- 
sponse to external fields. Nanoscale systems pose new 
challenges to the theory of strong correlations, since the 
latter are usually enhanced by spatial confinement. Vir- 
tually every future (nano) technology will use devices 
which interact with a time dependent (TD) environment. 
This increases the demand for ab initio methods to de- 
scribe realistic SCS acted upon by fast TD external fields. 

In the last decade, TDDFT [5] has emerged as an ef- 
fective ab initio treatment of TD phenomena [3J HJ [S]. 
DFT and TDDFT functionals, although related, are dif- 
ferent entities [6]: progress within TDDFT comes with 
progress with non equilibrium functionals. Constructing 
TDDFT functionals is an active area of research, with 
much work done, for example, in terms of the so-called 
Optimised Effective Method and extensions [3 [S]. A 
systematic route is given by a variational approach to 
Many Body Perturbation Theory (MBPT) [SJ, with a 
controlled improvement of the functionals [TD] . We also 
mention recent work |llj to include corrections to the 
ALDA [lj]. Current TDDFT functionals are quite suc- 
cessful for weakly interacting systems or in the linear 
response regime [3j HJ [5]. To date, no studies are avail- 
able of TDDFT applied to SCS (for DFT, see [UES]). 
At this early stage, model systems can be of aid, to pro- 
vide guidelines for ab initio approaches. An assessment 
of TDDFT for model SCS under TD fields is thus highly 
desirable. 

Here we study finite Hubbard chains in the presence of 
TD external fields, and use the results from exact time 
evolution to assess the potential of TDDFT for SCS. We 
also introduce an ALDA to v xc , based on an LDA-Bethc- 
Ansatz approach to the ground state of the inhomoge- 
neous ID Hubbard model [2]. Our main results are i) 
in the range of parameters we investigated, TDDFT is 
a practically viable route to describe SCS far away from 
equilibrium and in the TD regime; this is our central re- 



sult; ii) an exact analytic treatment for a two-site chain 
and an exact inequality for general ID chains are con- 
sistent with the numerical results; iii) for not too large 
external fields, the exact xc potential, v xc , obtained nu- 
merically by reverse engineering, is regular and well be- 
haved within the time span of our simulations. However, 
in some cases, v xc shows sharp structures in its temporal 
profile; iv) strong electron-electron interactions reduce 
memory effects; yet, non-adiabatic and non-local effects 
are in general necessary ingredients for a TDDFT of SCS. 
TDDFT time evolution for the many body problem. We 
study open-ended Hubbard chains, with Hamiltonian 

+fc(t)5>l, (1) 

(RR')a R <J 

with riRcr = a} Ra aR ai a =|,j and (R,R'} denoting 
nearest neighbour sites. The hopping parameter is V 
(V = —1), U is the interaction strength, and h(t) is the 
strength of a spin independent, local external field. U 
and h(t) are given in units of | . For simplicity, h(t) 
is localised at the leftmost site (R = 1), but we exam- 
ined other couplings, not discussed here. We consider 
L = 4 to 12 sites and N e = L or 3L/2 electrons (half- 
and three-quarter filling densities); we take spin up and 
down electrons equal in number; this holds during the 
time evolution, since H has no spin-flip terms. To evolve 
in time the exact many-body we use the Lanc- 

zos's algorithm [16 in the mid-point approximation (in 
all calculations the timestep A = 0.02|y| _1 ; numerical 
convergence was checked by halving A). By a fitting pro- 
cedure, we find an exact Kohn-Sham (KS) Hamiltonian, 

H KS = VY.(RR') a a+ ia a Rla +h{t)Y Ja n lcT + Y jRa v eff (R,t)n Ra , 

where v e ff(R,t) = vii(R,t)+v xc (R,t). Since we consider 
nonmagnetic regimes, v e ff(R,t) is spin-independent. We 
require the exact and the KS electron densities to be the 
same at each time and cluster site. In practice, v xc (R,t) 
is determined by minimising T.Rxj\{ n Ro)f S ~ { n R.<j)t\ 2 - 
Our procedure extends to TDDFT a method introduced 
long ago in DFT [17] . In [IS] , it was shown quite gener- 
ally how to map from TD densities to potentials. Such 
mapping was recently used for a He model atom |19j . 
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FIG. 1: (Color online) L = 4, N e = 4. Top left : the ex- 
ternal fields h(t) = h a ,b(t)- Mid and bottom rows: n(R,t) 
and v xc (R,t) for different (U,h) cases: black solid, orange 
solid (grey in b&w) and black dashed curves refer to sites 
R = 1,2,4, respectively (R=3 results can be obtained via 
J2 R Vxc{R,t) = 0,J2 R n{R,t) = 2). Top right: v xc /U at 
R — 1 for U — 3,hb (thin black curve), U = 15, h a , (or- 
ange/grey curve), U — 15, ht (thick black curve). All panels 
have the same time intervals/units. 

Half filled chains: Exact Results and TDDFT. Small clus- 
ters prevent local excitations from propagating away, 
introducing time oscillations in the density; as a way 
to mitigate size effects, we consider chains of different 
lengths. We start with a four-site chain at half-filling 
(Fig. [I]). In the initial, ground state, n{R,t = 0) = 
(g\nR a \g) — 0.5 at any site. We choose U = 3 and 15, as 
examples of two interaction regimes. The results for v xc 
are shown in the bottom panels: when v xc is reused in 
the KS equations, it reproduces n{R, t) = (nna-)t (middle 
panels) with an accuracy of 10~ 5 or better (this applies 
to all figures). Since v xc is defined up to an arbitrary 
site independent function C(t), we display the potentials 
differences, e.g. 5v xc (R,t) = v xc (R,t) — v™(£), where 
v xc(t) = (l/L)J2i v xc(Ri,t). This also applies to vr 
and v ef f( for the Hartree term, vf?(t) = UN e /2L). For 
simplicity, the prefix S will be omitted. Also, we find 
useful to rescale v xc ( and v e ff) by U when comparing 
results for different U's. In Fig. [T] we consider two ex- 
ternal perturbations, switched-on/off at a faster {h a ) or 
slower (hh) rate. For U — 15 and h = h a (t), the densities 
exhibit fast oscillations superimposed on a smoother, av- 
erage change. For U = 15 and hb(t), the oscillations are 
considerably suppressed, due to a more gradual change 
of the overlap between the initial, ground state and the 
excited ones during the onset of /i&(t). The degree of 
charge redistribution is determined by U : for example, 
for hb(t), the (small) charge imbalance at U = 15 is fully 
absorbed by the second, R — 2, site; at U — 3, all sites 
are involved. From a TDDFT perspective, this is a con- 
sequence of how v xc depends on U . In the bottom panels, 
we can see that v xc and n(R, t) behave rather similarly. 
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FIG. 2: (Color online) Half-filling results. Top left four pan- 
els: exact density (thick black) and v xc /U (thin black) curves 
at site R = 1 when (7 = 1 and L=6 to 12. The field h(£) is 
also shown (dashed curve). Right six panels: L = 8. Dashed 
(solid) curves refer to fti(t) = h(£) (h 2 (t) = 2t)(t)). Left bot- 
tom graphs: ground state LDOS at R = 1 for L = 8, with 
a site energy shift e R =i = 0,1,2. Darker (lighter) patterns 
refer to the hole (electron) LDOS. A Lorentzian broadening 
was introduced. 

On the other hand, results in the top right panel of Fig. 
[T]show that the range variation of the rescaled quantities, 
i.e. v xc /U, is comparable in all cases. Also, for large U, 
Vxc/U is very much in phase with the perturbation. For 
U = 3, out-of-phase effects are evident. For example, the 
largest oscillations in v xc /U occur when hb(t) has already 
returned to (this suggest that large U values tend to 
reduce memory effects). This is a rather generic behav- 
ior, that we noted also for other fields [15] . 

Larger chains at half filling. In Fig. [2] top left panel, we 
show results for L=6- to 12-site clusters, with [7=1 and 
h(t) = t)(t) the same for all L's. Chains with different L 
behave rather similarly the obvious differences being due 
to the fine details of the excited states. In Fig. [2] in the 
six panels on the right, we compare results for U = 1 and 
3 and two perturbations hi(t) = \)(t) and /i2(t) = 2f)(£), 
when L — 8. To get an idea of the strength of hi >2 {t), 
we can look at the equilibrium one-particle interacting 
local density of states, LDOS (Fig. [2] bottom left pan- 
els), when a static shift e_R-i = 1,2 is introduced. The 
shift corresponds to the maximum value h" l % x achieved 
by h\^{t) during the time evolution (the unshifted LDOS 
is also shown) and induces significant (larger for smaller 
U) spectral changes. The maxima = 1>2 are large 

enough to induce transitions from occupied to empty lev- 
els (cfr. with the energy gap in the unshifted LDOS). Re- 
sults for the TD density at U = 1 and U = 3 in the six 
right panels are consistent with the LDOS features and 
with results from Fig. [T] For example, charge variations 
are affected both by U and h(t): a larger U (a smaller 
field) induces a weaker response. Also, at U = 3, the TD 
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density imbalance is localised near the site R — 1; for 
U=l, it redistributes across all sites. Finally, in a broad 
parameter range, TDDFT reproduces the exact density. 
Adiabaticity vs locality and TDDFT. We now introduce 
an ALDA to TDDFT for the Hubbard model, and ap- 
ply it to a chain with L = 8 and N e — 12 (thus we 
also show how TDDFT performs away from half filling) . 
To disentangle adiabatic from locality effects in v xc , we 
use two approximations (Al and A2) for the density. In 
Al, we calculate at every timestep the ground state one- 
particle density of the instantaneous many body Hamil- 
tonian, Eq.Q (this implies no approximations based on 
local potentials). In A2, we introduce an ALDA to v xc ; 
our ALDA uses a a Local Density Approximation (LDA) 
to v xc for the ground state of the ID inhomogenous Hub- 
bard model [13], based on the Bethe-Ansatz (BA). We 



employ an analytical interpolation to v. 



BA-LDA 



(n) 



BALD A 



(n) = /j, [2cos(ttz//3) - 2cos(ttz/2) + Uz/2] (2) 



where n = n + + rt_, z = 1 — \n — 1|, (J, = sgn(n — 1) 
and f3 = f3(U), which is independent of n, is obtained 
from the BA solution at half filling [T3]. From Eq.([2|, 
v xc has a jump at n = 1, A xc = 4 cos(7r / /3(U)) + U. In 
our novel ALDA scheme, v x ALDA becomes a function 
of the instantaneous densities along the KS trajectories, 
= v x 3ALDA (n KS (R,t)). In Fig. [3J we compare ex- 
act and approximate results for slow and fast perturba- 
tions. All results are for site R = 1. We begin with the 
non-adiabatic case (panels in the top four rows) where 
U = 3,6,/ii(t) = E)(i) and h 2 (t) = 2f)(i), with f|(f) the 
same as m Fig. |2] An exact TDDFT description (black 
solid curves) is possible also away from half filling n > 1. 
As when n — 1, but to a lesser extent, larger U values 
reduce the changes in the density due to h(t). For /12(f) 
and U = 3, the exact v xc exhibits sharp resonances for 
t > 90, when /12(f) has returned to zero. At the same 
points, the density behaves smoothly. We observed such 
peaks for other kinds of perturbations and other param- 
eters values; their intensity increases at larger pertur- 
bation strengths. Such structures might be a challenge 
in constructing approximate potentials. Turning to Al, 
we note that a non-local but fully adiabatic description 
(dashed curves) gives correctly the average profile of the 
density. With the exact wavefunction replaced by the 
instantaneous (exact) ground state counterpart, there is 
no contribution from the excited states. This removes 
memory effects, and the density closely follows the tem- 
poral profile of h(t): for example, the oscillations in the 
exact densities around t — 90 are completely missed by 
Al. For A2, (orange curves, grey in b&w) the agreement 
with the exact densities is better, the maximum discrep- 
ancy being within a few percent (this level of discrep- 
ancy we also found for BALDA ground state densities). 
When /11.2(f) = 0), A2 reproduces many aspects of the 
exact results. However, a significant time-delay of certain 
traits suggests that memory effects are not being properly 



taken into account (this especially manifests for f > 90, 
when /11.2(f) = 0). The agreement between exact and A2 
densities looks better for /12(f) than for hi(t). To elabo- 
rate on this point, we compare exact and A2 results for 
v xc . For h\(t) and both U values, there is a small discrep- 
ancy between v xc and v A ^ . This accounts for part of the 
difference between the corresponding densities (the other 
part being due to vh)- For /12(f), when 2n KS (R,t) = 1, 
v A £ shows discontinuities (see Eq.([2])) which are absent 
in v xc (t). This suggests that, for h 2 (t), the agreement 
between A2 and exact results is somewhat accidental, 
whenever hit) drives n across the half-filling value. To 
corroborate this point, we show (Fig. [3j bottom row) re- 
sults for two slow perturbations hf(t) = f) a d.(f) (dashed 
curves) and /12(f) = 2f) a ^.(f) (solid curves), with fw.(f) 
a smoothcned version of f)(f) of Fig. [2] The exact (black 
curves) and Al densities are identical (Al densities are 
fully underneath the exact ones), i.e. the exact time evo- 
lution is fully adiabatic; A2 performs well (orange curves, 
grey in b&w) whenever n K (f) does not cross the half 
filling point. If the crossing occurs (U — 3, h?,(t)), we 
see noise-like features in n KS , due to the jump in v x ^. 
To summarize this section, TDDFT reproduces the ex- 
act density with a reasonable-looking v xc which can be, 
however, rather different from the ALDA one. And, in 
general, non-adiabatic and non-local effects (the jump is 
a non local feature) are both needed in v xc for a TDDFT 
of SCS. We finally note that for SCS with nearly filled 
bands the T-matrix approximation |20j . TMA, is quite 
successful [21] . As a mention of work in progress, a study 
of v xc in the TD TMA is under way. 

Remarks about v-representability. It was recently 
pointed out [52] that, for lattice models, there is an is- 
sue concerning the mapping of TD densities to poten- 
tials. The analysis in |22] was for a one-particle, two site 
system. To discuss the Kohn-Sham v-representability 
(KSVR) in a two-site system with N e = 1 electrons, we 
write, for any time f, ip K s(t) = e^(n 1 / 2 , e^(l - n) 1 / 2 ); 
one can show that \n\ < 2\V\(n(l — n)) 1 / 2 , a necessary 
and sufficient condition for KSVR. For the trial density 
in |22j . there are time intervals where such inequality 
is not fulfilled. Similar conclusions for L = 2 were in- 
dependently reached in [23], where the L > 2,N e = 1 
was also studied. In [23], the inequality for |n| was re- 
lated to the real vs complex nature of the effective po- 
tential. For L > 2, N e = 1 the condition for a real poten- 
tial is [23 j \S k \ < 2\V\ v /n k n k+1 \, where S k = J2i=i fk- 
However, for the present work, we need to consider the 
many-particle, interacting case. Starting with a Hub- 
bard dimer (HD), i.e. with Eq.([T|) for L = 2, we per- 
formed simulations for different pairs (U, h(t)) and veri- 
fied that \h\ < 2|V|(n(l - n)) 1 / 2 is always obeyed. This 
offers strong evidence of the KSVR of a HD. To com- 
plete a formal proof, one needs to show that such in- 
equality always holds for an HD. This is indeed the 
case [53]. The many-particle L > 2 case is consider- 
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FIG. 3: (Color online) L=8 at 3/4 filling. Top four rows: fast 
fields. Black solid, black dashed, and orange (grey in b&w) 
curves denote exact, Al and A2 results, respectively. The 
time unit is the same for all panels; panels in the same row 
share the same vertical scale. The curves in each row (i.e. 
Tii, Vxo/U) wee specified on the left, while the fields (i.e. hi or 
/12 ) are specified on the right. Bottom row: slow fields (note 
the horizontal axis). Solid (dashed) curves refer to hf (h^) 
fields. Black (orange/grey) curves denote exact (A2) results. 

ably more complicated, and here we limit our discus- 
sion to a simple but necessary condition for KSVR. For 
the KS system, the total density (per spin channel) at 
the A:-th site is n£. ot = J2\ n k' where A labels the KS 
one particle states. As a generalisation of the result 
in [23], we define Sg* = J2\ S k = Ei=i We § et 

\ s k 0t \ ^ J2\ 2 \ v \( n k n k+i) 1/2 and ' usin S the Schwarz in- 
equality, Sl ot < 2\V\{n t ° t n t °l 1 ) 1 / 2 . For the interacting 
many body system, we start with (n& s ) = i([H, Uks]) — 
iV([al +l s a k , s - a{ s a k - 1:S - h.c.}). We then get * = 
S i= i(^fes) = iV{[a\. +lta ak,s — h.c.]). By the same manip- 
ulations as in the HD,'|5f*| < 2\V\{n t ° t n t °l 1 ) 1 / 2 . Thus, 
for L > 2,N e > 1, the inequality holds for the KS and 
the interacting ID systems, which is consistent with the 
numerical results for L > 2. 

In conclusion, we provided a characterisation of TDDFT 
for strongly correlated systems. We compared exact vs. 
approximate results from the time evolution of model 
finite systems, in a broad range of model parameters. 
The exact v xc gave us insight into some of the prop- 
erties approximate xc functionals should satisfy. The 
v-represent ability problem was discussed, and an adia- 
batic approximation was introduced |25j . Our results 
illustrate the scope of TDDFT for non equilibrium phe- 
nomena in the presence of strong, time varying external 
fields in SCS, and encourage further investigations, some 
of which currently under way. We acknowledge many 
profitable discussions with C-O. Almbladh and U. von 
Barth. We also thank K. Capelle and K. Burke for use- 
ful conversations. This work was supported by the EU 
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